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Abstract 

We  have  previously  described  algorithms  for  a  battery  management  system  (BMS)  that  uses  Kalman  filtering  (KF)  techniques  to  estimate  such 
quantities  as:  cell  self-discharge  rate,  state-of-charge  (SOC),  nominal  capacity,  resistance,  and  others.  Since  the  dynamics  of  electrochemical  cells 
are  not  linear,  we  used  a  non-linear  extension  to  the  original  KF  called  the  extended  Kalman  filter  (EKF). 

We  were  able  to  achieve  very  good  estimates  of  SOC  and  other  states  and  parameters  using  EKF.  However,  some  applications  e.g.,  that  of 
the  battery-management- system  (BMS)  of  a  hybrid-electric-vehicle  (HEV)  can  require  even  more  accurate  estimates  than  these.  To  see  how  to 
improve  on  EKF,  we  must  examine  the  mathematical  foundation  of  that  algorithm  in  more  detail  than  we  presented  in  the  prior  work  to  discover  the 
assumptions  that  are  made  in  its  derivation.  Since  these  suppositions  are  not  met  exactly  in  BMS  application,  we  explore  an  alternative  non-linear 
Kalman  filtering  techniques  known  as  “sigma-point  Kalman  filtering”  (SPKF),  which  has  some  theoretical  advantages  that  manifest  themselves 
in  more  accurate  predictions.  The  computational  complexity  of  SPKF  is  of  the  same  order  as  EKF,  so  the  gains  are  made  at  little  or  no  additional 
cost. 

The  SPKF  method  as  applied  to  BMS  algorithms  is  presented  here  in  a  series  of  two  papers.  This  first  paper  is  devoted  primarily  to  deriving 
the  EKF  and  SPKF  algorithms  using  the  framework  of  sequential  probabilistic  inference.  This  is  done  to  show  that  the  two  algorithms,  which  at 
first  may  look  quite  different,  are  actually  very  similar  in  most  respects;  also,  we  discover  why  we  might  expect  the  SPKF  to  outperform  EKF  in 
non-linear  estimation  applications.  Results  are  presented  for  a  battery  pack  based  on  a  third-generation  prototype  LiPB  cell,  and  compared  with 
prior  results  using  EKF.  As  expected,  SPKF  outperforms  EKF,  both  in  its  estimate  of  SOC  and  in  its  estimate  of  the  error  bounds  thereof.  The 
second  paper  presents  some  more  advanced  algorithms  for  simultaneous  state  and  parameter  estimation,  and  gives  results  for  a  fourth-generation 
prototype  LiPB  cell. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

This  paper  applies  results  from  the  field  of  study  known  vari¬ 
ously  as  sequential  probabilistic  inference  or  optimal  estimation 
theory  to  advanced  algorithms  for  a  battery  management  sys- 
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tern  (BMS)  for  hybrid-electric  vehicle  (HEV)  application.  This 
BMS  is  able  to  estimate  battery  state-of-charge  (SOC),  instanta¬ 
neous  available  power,  and  parameters  indicative  of  the  battery 
state-of-health  (SOH)  such  as  an  increase  in  cell  resistance  (i.e., 
power  fade)  and  capacity  fade,  and  is  able  to  adapt  to  chang¬ 
ing  cell  characteristics  over  time  as  the  cells  in  the  battery  pack 
age.  The  algorithms  have  been  successfully  implemented  on  a 
lithium-ion  polymer  battery  (LiPB)  pack,  and  we  also  expect 
them  to  work  well  for  other  battery  chemistries. 

A  hybrid-electric  vehicle  is  one  with  both  an  internal- 
combustion  engine  and  an  electric  motor.  Both  may  be  cou¬ 
pled  directly  to  the  power  train — resulting  in  a  “parallel  hybrid” 
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configuration — where  the  motor  provides  boost  energy  to  sup¬ 
plement  the  engine,  and  acts  as  a  generator  when  coasting,  brak¬ 
ing,  or  when  the  engine  can  supply  extra  power  to  charge  the 
battery  pack.  Alternately,  the  engine  may  be  used  exclusively 
to  drive  a  generator  that  charges  the  battery  pack;  the  motor  is 
then  coupled  directly  to  the  power  train — resulting  in  a  “series 
hybrid”  configuration.  The  series  configuration  promises  greater 
potential  efficiency,  at  the  cost  of  a  larger  required  battery  pack. 
At  the  time  of  the  writing  of  this  paper,  the  only  HEVs  mar¬ 
keted  in  the  US  are  parallel  hybrid  systems  and  require  a  battery 
pack  of  fairly  modest  size.  Even  so,  and  because  demanding 
requirements  on  a  pack  of  limited  capacity  result  in  cell  electro¬ 
chemistries  that  are  often  far  from  equilibrium,  advanced  meth¬ 
ods  must  be  used  to  estimate  SOC,  SOH,  and  instantaneous 
power  in  order  to  safely,  efficiently  and  aggressively  exploit  the 
pack  capabilities. 

The  methods  we  use  to  estimate  these  numeric  quantities 
are  based  on  model-based-estimation  theory.  Two  separate  com¬ 
ponents  are  needed  to  implement  model-based  estimation:  (1) 
a  model  of  cell  dynamics,  and  (2)  an  algorithm  that  uses  that 
model.  The  general  approach  is  to  use  measured  model  inputs  to 
predict  measurable  model  outputs  using  present  values  of  model 
state  and  parameters.  Any  difference  between  the  model  out¬ 
put  and  the  measured  system  output  can  be  attributed  to  errors 
in  the  measurements,  states,  parameters,  or  the  model  frame¬ 
work  itself.  An  algorithm  uses  this  difference  signal  to  intelli¬ 
gently  update  its  estimate  of  the  model  state  and/or  parameters. 
In  this  work,  the  model  of  cell  dynamics  is  the  “enhanced  self- 
correcting”  (ESC)  model,  introduced  elsewhere  [3,5]  and  briefly 
reviewed  in  Section  7.2  here.  However,  the  model  is  not  the  fo¬ 
cus  of  this  work,  but  rather  the  algorithms  that  use  that  model, 
which  comprise  the  family  of  Kalman  filters.2 

Kalman  filters  are  an  intelligent — and  sometimes  optimal — 
means  for  estimating  the  present  value  of  the  time-varying 
“state”  of  a  dynamic  system.  By  modeling  our  battery  system  to 
include  the  wanted  unknown  quantities  in  its  state  description, 
we  may  use  a  Kalman  filter  to  estimate  their  values.  An  addi¬ 
tional  benefit  of  the  Kalman  filter  is  that  it  automatically  provides 
dynamic  error-bounds  on  these  estimates  as  well.  We  exploit  this 
fact  to  give  aggressive  performance  from  our  battery  pack,  with¬ 
out  fear  of  causing  damage  by  overcharge  or  overdischarge. 

We  have  previously  reported  work  using  extended  Kalman 
filters  (EKF)  to  solve  the  BMS  algorithm  requirements  [1-6]. 
We  have  since  explored  a  different  form  of  Kalman  filtering 
called  sigma-point  Kalman  filters  (SPKF),  and  have  found  them 
to  have  several  important  advantages  to  be  outlined  herein. 

We  present  the  base-line  SPKF  and  some  variants,  along  with 
testing  results  and  analysis,  in  a  two-part  series  of  papers.  This 
first  paper  primarily  derives  the  equations  that  govern  optimal 
Gaussian  sequential  probabilistic  inference,  and  then  shows  how 
these  equations  apply  to  standard  Kalman  filtering,  extended 
Kalman  filtering,  and  sigma-point  Kalman  filtering.  The  SPFK 


2  These  algorithms  may  be  applied  to  any  cell  model  in  the  correct  format,  and 
hence  may  be  used  to  estimate  the  state  and  parameters  of  cells  with  differing 
electrochemistries  and  physical  configurations. 


algorithm  is  exercised  on  the  same  data  as  presented  in  [5]  to 
compare  with  EKF  on  a  third-generation  prototype  LiPB  HEV 
cell.  In  all  cases,  the  SPKF  outperforms  the  EKF,  both  in  its 
ability  to  estimate  SOC,  and  in  its  estimates  of  the  error  bounds 
thereof. 

The  second  paper  in  this  series  [7]  explores  some  more  ad¬ 
vanced  algorithms:  square-root  sigma-point  Kalman  filtering, 
parameter  estimation,  and  simultaneous  state  and  parameter  es¬ 
timation  using  the  joint  and  dual  methods.  Since  it  has  already 
been  established  in  this  first  paper  that  SPKF  outperforms  EKF, 
we  no  longer  use  data  from  the  (now  obsolete)  third-generation 
cells,  but  present  testing  data  from  fourth-generation  FiPB  high- 
power  HEV  cells. 

While  there  is  necessarily  some  redundancy  in  presenting 
this  work  in  a  two-paper  series,  we  felt  that  this  approach  had 
certain  advantages  to  a  single  monolithic  manuscript.  First,  each 
part  is  self-contained,  with  the  first  paper  primarily  serving  to 
show  the  theoretical  benefits  of  SPKF  and  to  demonstrate  these 
with  examples  and  the  second  paper  primarily  serving  to  show 
how  the  SPKF  might  be  used  in  actual  application;  secondly, 
the  results  of  the  first  paper  are  found  for  the  third-generation 
cell,  and  are  directly  comparable  to  previously  reported  results, 
while  the  results  of  the  second  paper  are  found  for  a  fourth- 
generation  cell  only — since  we  have  already  established  that  the 
SPKF  outperforms  EKF,  there  is  no  need  to  rehash  old  data  with 
the  SPKF  and  burden  this  work  with  an  excess  of  figures  and 
tables — rather,  we  can  lay  the  foundation  for  future  work  that 
will  be  reported  for  these  cells. 

2.  Sequential  probabilistic  inference 

Very  generally,  any  causal  dynamic  system  (e.g.,  a  battery 
cell)  generates  its  outputs  as  some  function  of  its  past  and  present 
inputs.  Often,  we  can  define  a  state  vector  for  the  system  whose 
values  together  summarize  the  effect  of  all  past  inputs.  Present 
system  output  is  a  function  of  present  input  and  present  state 
only;  past  input  values  need  not  be  stored.  The  system’s  pa¬ 
rameter  vector  comprises  all  quasi-static  numeric  quantities  that 
describe  how  the  system  state  evolves  and  how  the  system  out¬ 
put  may  be  computed.  The  state- vector  quantities  change  on  a 
relatively  rapid  time  scale,  and  the  parameter-vector  quantities 
change  on  a  relatively  long  time  scale  (or,  not  at  all). 

For  some  applications,  we  desire  to  estimate  the  state-  or 
parameter- vector  quantities  in  real  time,  as  the  system  operates. 
Probabilistic  inference  and  optimal  estimation  theory  are  names 
given  to  the  field  of  study  that  concerns  estimating  these  hidden 
variables  in  an  optimal  and  consistent  fashion,  given  noisy  or 
incomplete  observations.  For  example,  in  this  paper  we  will  be 
concerned  with  estimating  the  state  of  an  electrochemical  cell 
and  the  parameters  for  a  mathematical  model  describing  the  cell 
dynamics.  SOC  is  one  element  of  the  state  vector  of  particular 
concern,  and  factors  such  as  resistance  increase  and  power  fade 
are  elements  from  the  parameter  vector  of  concern  to  help  esti¬ 
mate  SOH.  Observations  are  available  to  us  at  sampling  points 
and  include:  cell  current  /&,  cell  terminal  voltage  yk,  and  cell 
temperature  7^,  where  the  subscript  k  indicates  that  the  mea¬ 
surement  is  taken  at  the  kth  sampling  point. 
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In  the  following,  we  will  assume  that  the  electrochemical 
cell  under  consideration  may  be  modeled  using  a  discrete-time 
state- space  model  of  the  form 

Xk  =  f(xk- 1,  Uk- 1,  Wk-uk  -  1)  (1) 

yk  =  h(xk,  uk,vk,  k).  (2) 

Here,  xk  G  W1  is  the  system  state  vector  at  time  index  k, 
and  Eq.  (1)  is  called  the  “state  equation”  or  “process  equation”. 
The  state  equation  captures  the  evolving  system  dynamics. 
System  stability,  dynamic  controllability  and  sensitivity  to 
disturbance  may  all  be  determined  from  this  equation.  The 
known/deterministic  input  to  the  system  is  ^  G  F,  and 
wk  G  M77  is  stochastic  “process  noise”  or  “disturbance”  that 
models  some  unmeasured  input  which  affects  the  state  of  the 
system.  The  output  of  the  system  is  yk  G  Mm,  computed  by 
the  “output  equation”  (2)  as  a  function  of  the  states,  input,  and 
Vk  G  Mm,  which  models  “sensor  noise”  that  affects  the  measure¬ 
ment  of  the  system  output  in  a  memoryless  way,  but  does  not 
affect  the  system  state.  f(xk- 1,  uk~  i,  wk~  i,  k  —  1)  is  a  (possi¬ 
bly  non-linear)  state  transition  function  and  g(xk,  uk,  vk,  k)  is 
a  (possibly  non-linear)  measurement  function. 

With  this  model  structure,  the  evolution  of  unobserved  states 
and  observed  measurements  may  be  visualized  as  shown  in 
Fig.  1.  The  conditional  probability  p(xk \xk-i)  indicates  that  the 
new  state  is  a  function  of  not  only  the  deterministic  input  uk- 1, 
but  also  the  stochastic  input  wk~  \ ,  so  that  the  unobserved  state 
variables  do  not  form  a  deterministic  sequence.  Similarly,  the 
conditional  probability  density  function  p(yk\xk)  indicates  that 
the  observed  output  is  not  a  deterministic  function  of  the  state, 
due  to  the  stochastic  input  vk. 

The  goal  of  probabilistic  inference  is  to  create  an  estimate  of 
the  system  state  given  all  observations  Yk  =  {yo,  yi ,  ■  •  • ,  yk}-  A 
frequently  used  estimator  is  the  conditional  mean 


xk  =  E[**|Y*]  = 


xkp(xk\Yk)dxk, 


where  RXk  is  the  set  comprising  the  range  of  possible  xk ,  and  E[-] 
is  the  statistical  expectation  operator.  The  optimal  solution  to  this 
problem  computes  the  posterior  probability  density  p(xk\Yk)  re¬ 
cursively  with  two  steps  per  iteration  [8] .  The  first  step  computes 
probabilities  for  predicting  xk  given  all  past  observations 


p(xk\Yk-i) 


p(xk  \xk-i)  p(xk- 1 1 ' Yk- 1 )  dxk- 1 


and  the  second  step  updates  the  prediction  via 


p(xk\Yk ) 


p(yk\xk)p(xk\Yk-i) 

p(yk\Yk-i) 


Fig.  1.  Sequential  probabilistic  inference. 


which  is  a  simple  application  of  Bayes’  rule  and  the  assumption 
that  the  present  observation  yk  is  conditionally  independent  of 
previous  measurements  given  the  present  state  xk.  The  relevant 
probabilities  may  be  computed  as 


p(yk\^k-i) 


p(yk\xk)p(xk\Yk~i)  dxk 


p{Xk\Xk- 1) 


Y  p(u» 

{w:xk=f(xk-i,uk-i,w,k—l)} 


p(yk\xk)=  Y  pW- 

{v:yk=h(xk,Uk,v,k)} 

Although  this  is  the  optimal  solution,  finding  a  formula 
solving  the  multi-dimensional  integrals  in  a  closed-form  is  in¬ 
tractable  for  most  real-world  systems.  For  applications  that  jus¬ 
tify  the  computational  expense,  the  integrals  may  be  closely  ap¬ 
proximated  using  Monte  Carlo  methods  such  as  particle  filters 
[9-1 2,8] .  However,  for  battery  packs  containing  tens  or  hundreds 
of  cells,  and  hence  perhaps  hundreds  or  thousands  of  states,  the 
present  economics  do  not  make  this  option  feasible. 

A  simplified  solution  to  these  equations  may  be  obtained  if 
we  are  willing  to  make  the  assumption  that  all  probability  den¬ 
sities  are  Gaussian — this  is  the  basis  of  the  original  Kalman 
filter,  the  extended  Kalman  filter,  and  the  sigma-point  Kalman 
filters  to  be  discussed.  Then,  rather  than  having  to  propagate  the 
entire  density  function  through  time,  we  need  only  to  evaluate 
the  conditional  mean  and  covariance  of  the  state  (and  parame¬ 
ters,  perhaps)  once  each  sampling  interval.  It  can  be  shown  (cf. 


Appendix  A)  that  the  recursion  becomes: 

xf  =  xj  +  Lk(yk  -  %)  =  xk  +  Lkyk  (3) 

Y,k  =  —  (4) 

where  the  superscript  T  is  the  matrix/vector  transpose  operator, 
and 

x\  =  E[xjt|Yjt]  (5) 

xk  =E[xk\Yk-i]  (6) 

fc  =  E[yjt|Y*_i]  (7) 

T,yk  =  E[(xk  -  xk)(xk  -  xkf]  =  E[(x^)(x^)t]  (8) 

Y,~xp  =  E[(xk  -  xf(xk  -  xfT]  =  E[(x^)(x^)t]  (9) 

=  E[(ji  -  %)(yk  -  % )T]  =  E[(5x)(5’(0T]  (10) 

Lk  =  E  [(xk  -  x~k)(yk  -  »)T]Sd  =  Sr.  (1 1) 


While  this  is  a  linear  recursion,  we  have  not  directly  assumed  that 
the  system  model  is  linear.  In  the  notation  we  use,  the  decoration 
“circumflex”  indicates  an  estimated  quantity  (e.g.,  x  indicates 
an  estimate  of  the  true  quantity  x).  A  superscript  “— ”  indicates 
an  a  priori  estimate  (i.e.,  a  prediction  of  a  quantity’s  present 
value  based  on  past  data)  and  a  superscript  “+”  indicates  an 
a  posteriori  estimate  (e.g.,  x ^  is  the  estimate  of  true  quantity 
v  at  time  index  k  based  on  all  measurements  taken  up  to  and 
including  time  k).  The  decoration  “tilde”  indicates  the  error  of 
an  estimated  quantity.  The  symbol  Yxy  =  E[xyT]  indicates  the 
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auto-  or  cross-correlation  of  the  variables  in  its  subscript.  (Note 
that  often  these  variables  are  zero-mean,  so  the  correlations  are 
identical  to  covariances.)  Also,  for  brevity  of  notation,  we  often 
use  to  indicate  the  same  quantity  as  Hxx. 

In  the  following  sections  we  will  apply  Eqs.  (3)— (1 1),  and 
approximations  thereof,  with  different  sets  of  assumptions  to 
derive  the  Kalman  filter,  the  extended  Kalman  filter,  and  sigma- 
point  Kalman  filters.  All  members  of  this  family  of  filters  will 
comply  with  a  structured  sequence  of  six  steps  per  iteration,  as 
outlined  here. 

General  step  1:  State  estimate  time  update.  Each  measure¬ 
ment  interval,  the  first  step  is  to  compute  an  updated  prediction 
of  the  present  value  of  Xk,  based  on  a  priori  information  and  the 
system  model.  This  is  done  using  Eqs.  (1)  and  (6)  as 

xk  =  E[jc*|Y*_i]  =  E[f(xk-i,Uk-i,  Wk~\ ,k  -  l)|Yfc_i], 

General  step  2:  Error  covariance  time  update.  The  second  step  is 
to  determine  the  predicted  state-estimate  error  covariance  matrix 
'Z%k  based  on  a  priori  information  and  the  system  model. 

We  compute  5^T k  =  E[(v^)(v^)T]  using  Eq.  (8),  knowing 
that  x^  =  Xk  —  x^ . 

General  step  3:  Estimate  system  output  yk.  The  third  step  is  to 
estimate  the  system’s  output  using  present  a  priori  information 
and  Eqs.  (2)  and  (7) 

%  =  E[yjfc|Y*_i]  =  E[h(xk,  uk,  vk,  fc)|Y*_i.] 


General  step  h-:  rssumaior  gain  matrix  i^k-  me  louru 

compute  the  estimator  gain  matrix  Lk  by  evaluating 
1 

knx  k' 


is  to 

^ xy,k  ^ y,k  ’ 


General  step  5:  State  estimate  measurement  update.  The  fifth 
step  is  to  compute  the  a  posteriori  state  estimate  by  updating 
the  a  priori  estimate  using  the  estimator  gain  and  the  output 
prediction  error  yk  —  %  using  (3).  There  is  no  variation  in  this 
step  in  the  different  Kalman  filter  methods;  implementational 
differences  between  Kalman  approaches  do  manifest  themselves 
in  all  other  steps,  however. 


Table  1 

Summary  of  the  general  sequential  probabilistic  inference  solution 

General  state-space  model 

Xk  =  f{xk-\,Uk-\,  Wk-i ,k  -  1) 
yk  =  h(xk,  uk,  vk,  k ), 

where  wk  and  vk  are  independent,  Gaussian  noise  processes  of  covariance 
matrices  Ew  and  E„  respectively 

Definitions:  let 

xf  -  xk  -xf,  yk  =  yk-  % 

Initialization:  for  k  =  0,  set 

Xq  =  E[*0] 

=  E[(*0  -  *o  )Oo  -  ^o)T] 

Computation:  for  k  =  1,2,...  compute 

State  estimate  time  update:  xf  =  E[/(x&_i,  uk- 1,  wk- 1,  k  —  1)|  Yjt— i] 
Error  covariance  time  update:  ET k  =  E[(v^  )(v^)T] 

Output  estimate:  %  =  E [h(xk,  uk,  vk,  &)|Yfc_i] 

Estimator  gain  matrix:  Lk  =  E[(v^)(yyt)T](IE[(y^)(y^)T])~1 
State  estimate  measurement  update:  x ^  =  xf  +  Lk(yk  —  yk) 

Error  covariance  measurement  update:  E fk  =  E2^  —  Lk'Ey^Lj 


General  step  6:  Error  covariance  measurement  update.  The 
final  step  computes  the  a  posteriori  error  covariance  matrix  using 
Eq.  (4).  The  estimator  output  comprises  the  state  estimate  x ^  and 
error  covariance  estimate  The  estimator  then  waits  until 
the  next  sample  interval,  updates  k,  and  proceeds  to  step  1 . 

The  general  sequential  probabilistic  inference  solution  is 
summarized  in  Table  1.  The  following  sections  describe  spe¬ 
cific  applications  of  this  general  framework.  First,  we  derive  the 
Kalman  filter,  and  then  the  extended  Kalman  filter.  This  is  done 
to  give  confidence  that  the  sequential  probabilistic  inference  ap¬ 
proach  does  indeed  result  in  the  equations  for  these  standard 
algorithms,  to  show  the  assumptions  made  when  performing  the 
derivations,  and  to  give  insight  into  how  we  might  improve  es¬ 
timation  for  non-linear  systems.  The  reader  may  wish  to  skip 
ahead  to  Section  5  on  the  first  reading  of  this  paper,  and  then 
return  to  Sections  3  and  4  for  a  greater  depth  of  understanding. 

3.  Optimal  application  to  linear  systems:  the  Kalman 
filter 


In  this  section,  we  take  the  general  results  from  Section  2 
and  apply  them  to  the  specific  case  where  the  system  dynam¬ 
ics  are  linear.  Linear  systems  have  the  desirable  property  that 
all  probability  distributions  do  in  fact  remain  Gaussian  if  the 
stochastic  inputs  are  Gaussian,  so  the  assumptions  made  in  de¬ 
riving  the  filter  steps  hold  exactly.  If  the  system  dynamics  are 
linear,  then  the  Kalman  filter  (first  presented  in  [  1 3 , 1 4] )  is  the  op¬ 
timal  minimum-mean- squared- error  and  maximum-likelihood 
estimator.  The  format  of  this  section  is  to  first  introduce  the 
form  of  a  linear  state-space  model,  and  then  to  apply  the  six 
steps  from  Section  2  to  this  form  to  derive  the  linear  Kalman 
filter  equations. 

The  linear  Kalman  filter  assumes  that  the  system  being  mod¬ 
eled  can  be  represented  in  the  “state-space”  form 

Xk  =  Ak-iXk-i  +  Bk-iup-i  +  Wk- 1 


yk  —  GkXk  T  DkUk  T  Vk. 


The  matrices  Ak  e  Wlxn,  Bk  e  Wlxp,  Ck  E  Wnxn  and  Dk  E 
Wnxp  describe  the  dynamics  of  the  system,  and  are  possibly 
time  varying.  Also,  both  Wk  and  Vk  are  assumed  to  be  mutually 
uncorrelated  white  Gaussian  random  processes,  with  zero  mean 
and  covariance  matrices  with  known  value: 


n  =  k\ 
0,  n  ^  k. 


'Ey,  n  =  k\ 
0,  n  ^  k. 


The  assumptions  on  the  noise  processes  Wk  and  Vk  and  on  the 
linearity  of  system  dynamics  are  rarely  (never)  met  in  practice, 
but  the  consensus  of  the  literature  and  practice  is  that  the  method 
still  works  very  well. 

KF  step  1:  State  estimate  time  update.  Here,  we  compute 


xk  =  E[Ak-\Xk-i  +  Bk-iUk-i  +  Wk- i|Yfc_i] 

—  Ak—\Xji_Y  T  Bk—  i  Uk—  i , 

by  the  linearity  of  expectation,  noting  that  Wk-i  is  zero-mean. 

KF  step  2:  Error  covariance  time  update.  First,  we  note 
that  the  estimation  error  may  be  found  by  comparing  Xk  = 
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Table  2 

Summary  of  the  linear  Kalman  filter  from  reference  [15] 

Linear  state-space  model 
xk  =  Ak- \Xk-i  +  Bk-\Uk-\  +  wk- 1 
yk  =  Ckxk  +  Dkuk  +  vk, 

where  wk  and  vk  are  independent,  zero-mean,  Gaussian  noise  processes 
of  covariance  matrices  Ew  and  E^,  respectively 

Initialization:  for  k  =  0,  set 
Xq  =  E[v0] 

EtQ  =  E[(*0  ~  *o  )Oo  “  ^o)T] 

Computation:  for  k  =  1,2,...  compute 

State  estimate  time  update:  xf  =  Ak- +  Bk-\uk-i 
Error  covariance  time  update:  E~  k  =  Ak-\ ET^_1  Aj,_[  +  E^ 

Output  estimate:  %  =  Ckxk  +  Dkuk 

Estimator  gain  matrix:  Lk  =  Ef  kCj [Ck  Ef  kCj  +  EJ-1 

State  estimate  measurement  update:  xk  =  xf  +  Lk(yk  —  yk) 

Error  covariance  measurement  update:  Efk  =  (7  —  LkCk) Ef  k 


Ak-iXk-i  +  Bk-iUk-i  +  Wk- 1  with  xk  as  computed  in  step  1. 
We  find  that  x k  =  Ak~\xk_x  +  Wk-\ .  Therefore, 

Er  =  E  l(xk)(~xkf] 

=  E[(A^_ixLi  +  u>k-i)(Ak-ixf  i  +  u^-i)7] 

=  Ak-iYfk_^Ak_i  + 

The  cross  terms  drop  out  of  the  final  result  since  the  white  state 
noise  Wk-i  is  not  correlated  with  the  state  at  time  k  —  1. 

KF  step  3:  Estimate  system  output.  We  estimate  the  system 
output  as 

yk  =  EfC^x^  +  DkUk  T  =  CkXk  T  E)kUk, 

since  Vk  is  zero-mean. 

KF  step  4:  Estimator  ( Kalman )  gain  matrix.  To  compute  Lk , 
we  first  need  to  compute  several  covariance  matrices.  Since  we 
know  that  yk  =  CkXk  +  F>kUk  +  Vk,  it  follows  that  %  =  + 

Vk  and 


4.  An  approximation  for  non-linear  systems:  the 
extended  Kalman  filter 


The  extended  Kalman  filter  is  one  approach  to  general¬ 
ize  the  KF  results  to  non-linear  systems.  At  each  point  in 
time,  steps  2  and  4  linearize  the  non-linear  state  and  output 
equations  around  their  present  operating  point  using  Taylor- 
series  expansions.  Steps  1  and  3  approximate  the  a  priori 
state  estimate  and  output  estimate  using  previously  computed 
terms. 

EKF  step  1:  State  estimate  time  update.  The  state  prediction 
step  is  approximated  as 


xk  =  E[/(x*_i,  w*_i,  wk-\ ,k  -  l)|Y*_i] 

%  f(x k_v  uk- 1,  wk-\,k-  1), 

where  Wk-i  =  E[u>fc_i].  That  is,  we  approximate  the  expected 
value  of  the  new  state  by  assuming  that  it  is  reasonable  to  sim¬ 
ply  propagate  xk_x  and  Wk-i  through  the  state  equation.  Often, 
Wk- 1  =  0. 

EKF  step  2:  Error  covariance  time  update.  The  covariance 
prediction  step  is  accomplished  by  first  making  an  approxima¬ 
tion  for  x f  . 

X^  =  Xk  -  Xk  =  f(xk-l,  Uk-1,  wk- I,k  -  1) 

-  f(Xk-VUk- 1,  Wk- L  k  -  1). 


The  second  term  is  expanded  as  a  Taylor  series  around 
the  prior  operating  “point”  which  is  the  set  of  values 

{xk-i,  uk-i,  Wk- 1,  k-  1} 


xk  ~ 


+ 


f(Xk- 1,  Uk-1,  Wk-I,k  -  1) 
df(xk-i,Uk-i,wk-uk  -  1) 


dxk-i 


(xfi  ~xk-i) 


- v - 

y\ 

Defined  as  Ak-i 


xk-i=4~\ 

✓ 


^y,k  —  Ck^x,k^k  + 

Again,  the  cross  terms  are  zero  since  Vk  is  uncorrelated  with  the 
estimate  x k  .  Similarly, 


+ 


df(xk-i,  Uk-i,  wk~\,  k  -  1) 
dwk-i 


(Wk-l  ~  wk-i). 


- V - 

Defined  as  Bk-i 


Wk-l=Wk-l 

✓ 


EVk  yj]  =  Mxk  (Ckxk  +  U£)T]  = 

Combining, 

Lk  =  ^kcTk[ck^kcJ  + 1:,]-1 . 

KF  step  6:  Error  covariance  measurement  update.  Finally, 
we  update  the  error  covariance  matrix. 

s %k  =  s %k  -  LkMykylvl 

=  -  Lk^[ykykWl%yk])~Ti&[xkyk]T 

~  ^x,k  ~  LkCk^x,k  =  (I  ~  LkCk)^x,k' 

The  final  Kalman  filter  for  linear  systems  is  summarized  in  Table 
2.  These  are  the  same  equations  as  presented  (without  derivation) 
in  reference  [4] . 


_  Z\  z\ 

This  gives  xk  ~  Ak- lXp-i  +  Bk-iWk-i-  Substituting  this  to 
find  the  predicted  covariance: 


^x,k  ~  )(*fc  )T]  ~  Ak-l^k-lAk-l  +  Bk-l'EwBk_1. 


+ 


,T 


EKF  step  3:  Output  estimate.  The  system  output  is  estimated 
to  be 


yk  =  E [h(xk,  uk,  vk,  fc)|Y*_i]  «  h(xk  ,  uk,  Vk,  k), 

where  Vk  =  E [vk].  That  is,  it  is  assumed  that  propagating  xk 
and  the  mean  sensor  noise  is  the  best  approximation  to  estimating 
the  output. 

EKF  step  4:  Estimator  gain  matrix.  The  output  prediction 
error  may  then  be  approximated 

%  =  yk  -  yk  =  h(xk,  Uk,  Vk,  k)  -  h(xk,  Uk,  Vk,  k) 
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Table  3 

Summary  of  the  non-linear  extended  Kalman  filter  from  reference  [16] 


Non-linear  state-space  model 
Xk  =  f(Xk- 1,  Wjfc-l,  Wk- i,k  -  1) 
yk  =  h(xk,  uk,  vk,  k ), 

where  wk  and  vk  are  independent,  Gaussian  noise  processes  of  covariance 
matrices  Eg,  and  E y,  respectively 


Definitions 

2,  _  df(xk,uk,wk,k) 

d  xk  Xk=xt 

q  _  dh(xk,uk,vk,k) 

d xk  xk  =xk 

Initialization:  for  k  =  0,  set 


Bk  = 

/V 

Dk  = 


*0  =  E[x0] 

S+0  =  E[(x0  -  xj)(*0  -  x+)T] 


df(xk,uk,wk,k ) 
dwk 

dh(xk,uk,vk,k) 

dvk 


wk=wk 
Vk=v  k 


Computation:  For  k  =  1,2,...  compute 

State  estimate  time  update:  xf  =  /(v^~_  j ,  i ,  k  —  1) 

Error  covariance  time  update:  E~  /.=A/:_iET^_1A^_1+JS^_iEu,JS^_1 
Output  estimate:  %  =  h(x7 ,  uk,  vk,  k ) 

_  A  rp  A  _  A  rp  A  _  /V  rp  1 

Estimator  gain  matrix:  Lk  =  E~  kCk  [C&E~  k Ck  +  Dk  EkDk  ] 

State  estimate  measurement  update:  xk  =  xf  +  Lk(yk  —  yk) 

I  ZV  _ 

Error  covariance  measurement  update:  ET^  =  (7  —  EkCk) E~  k 


using  again  a  Taylor-series  expansion  on  the  second  term. 

3 h(xk,  uk,  vk,  k) 


%  ^  h(xk,uk,  vk,k)  + 


3 xk 


(xk  -xk) 


xk=x. 


- V - 

-A, 

Defined  as  Ck 


+ 


3 Kxk,  uk,  k) 

9  vk 


(Vk  -  vk). 


V* 

/V 

Defined  as 


From  this,  we  can  compute  such  necessary  quantities  as 

«  QEr.cJ  +  hfcEshJ, 

^ly.k  ^  Wxk)(CkXk  +  bkVk)T]  =  ^ikCj- 

These  terms  may  be  combined  to  get  the  Kalman  gain 

Lk  =  ^kCJk[CkY7kCJk  +  Dk^vDTkrl. 

EKF  step  6:  Error  covariance  measurement  update. 

Finally,  the  updated  covariance  is  computed  as 

^ x,k  ~  Lk'EyfLj  =  E~  ^  -  Lk'Ey,k('Ey,k)  T(^xy,k^ 

=  ^Xf  ~  k  =  (/  -  LkCk)^y.  k. 


This  completes  the  derivation  of  the  extended  Kalman  filter. 
It  is  summarized  in  Table  3.  These  are  the  same  equations  as 
presented  (without  derivation)  in  reference  [4]. 


5.  Problems  with  the  EKF 


calculation  of  the  output  random  variable  mean,  the  other  con¬ 
cerns  the  output  random  variable  covariance. 

First,  we  note  that  EKF  step  1  attempts  to  determine  an  output 
random- variable  mean  from  the  state-transition  function  /(•) 
assuming  that  the  input  state  is  a  Gaussian  random  variable. 
EKF  step  3  makes  a  similar  calculation  for  the  output  function 
h(-).  The  EKF  makes  the  simplification 

E[fn(x)]  ~  fn(E[x]), 

which  is  not  true  in  general,  and  not  necessarily  even  close  to 
true  (depending  on  “how  non-linear”  the  function  fn(-)  is).  The 
SPKF  to  be  described  will  make  an  improved  approximation  to 
the  means  in  steps  1  and  3. 

Secondly,  in  EKF  steps  2  and  4,  a  Taylor-series  expansion  is 
performed  as  part  of  a  calculation  designed  to  find  the  output- 
variable  covariance.  Non-linear  terms  are  dropped  from  the  ex¬ 
pansion,  resulting  in  a  loss  of  accuracy.  The  SPKF  uses  a  dif¬ 
ferent  method  to  compute  covariances  and  will  improve  these 
estimates  as  well. 

To  give  a  simple  one-dimensional  example  illustrating  these 
two  effects,  consider  Fig.  2.  The  non-linear  function  is  drawn, 
and  the  input  random- variable  PDF  is  shown  on  the  horizontal 
axis,  with  mean  1.05.  The  straight  dotted  line  is  the  linearized 
approximation  used  by  the  EKF  to  find  the  output  mean  and 
covariance.  The  output  approximate  PDF  estimated  by  EKF  is 
drawn  as  a  dotted  line  on  the  vertical  axis,  where  a  Gaussian  PDF 
with  the  same  mean  and  variance  of  the  true  data  is  shown  as 
a  solid  PDF  on  the  same  axis.  We  notice  significant  differences 
between  the  means  and  covariances,  indicating  that  EKF  is  not 
producing  an  accurate  estimate  of  either  one. 

For  a  two-dimensional  example,  consider  Fig.  3.  Frame  (a) 
shows  a  cloud  of  Gaussian-distributed  random  points  used  as 
input  to  this  function,  and  frame  (b)  shows  the  transformed  set  of 
output  points.  The  actual  95%  confidence  interval  (indicative  of 
a  contour  of  the  Gaussian  PDF  describing  the  output  covariance 
and  mean)  is  shown  as  a  solid  ellipse  with  the  output  mean 


EKF  versus  SPKF  mean  and  variance  approximation 


Function  input  value 


The  extended  Kalman  filter  is  probably  the  best  known  and 
most  widely  used  non-linear  Kalman  filter.  However,  it  has  a 
number  of  flaws  that  can  be  improved  upon  fairly  easily  to  im¬ 
prove  state  estimation.  These  flaws  reside  in  two  assumptions 
made  in  order  to  propagate  a  Gaussian  random  state  vector  a 
through  some  non-linear  function:  one  assumption  concerns  the 


Fig.  2.  EKF  vs.  SPKF  mean  and  variance  approximation.  The  solid  line  is  the 
non-linear  function;  the  dotted  straight  line  is  the  linear  approximation  used 
by  EKF  in  the  neighborhood  of  input  variable  equal  to  1.05.  The  solid  squares 
on  the  axes  are  the  input  and  output  sigma  points.  The  solid-line  PDFs  are  the 
Gaussian  PDFs  matching  the  input  and  matching  the  output  mean  and  variance; 
the  dashed-line  PDF  is  the  output  PDF  predicted  by  the  sigma-point  method, 
and  the  dotted  PDF  is  the  output  PDF  predicted  by  the  EKF  method. 
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Input  data  scatter 


Fig.  3.  Two-dimensional  Gaussian  random  data  (gray  points  in  (a))  processed  by 
a  non-linear  function  to  become  the  gray  points  in  (b).  The  solid  ellipses  are  the 
true  95%  confidence  bounds,  the  black  squares  are  the  input  and  output  sigma 
points;  the  dashed  ellipse  is  the  95%  confidence  bound  produced  by  the  sigma- 
point  method,  and  the  dotted  ellipse  is  the  95%  confidence  bound  produced  by 
the  EKF  method. 

being  at  the  center  of  the  ellipse.  The  dotted  ellipse  shows  the 
covariance  predicted  by  EKF,  with  the  EKF  mean  being  at  the 
center  of  that  ellipse.  Again,  EKF  is  very  far  from  the  truth. 

In  both  examples,  SPKF  greatly  outperforms  EKF.  We  dis¬ 
cuss  why  in  the  next  section. 

6.  Sigma-point  Kalman  filters 

We  have  seen  that  the  EKF  approach  to  generalizing  the  KF 
to  non-linear  systems  is  to  linearize  the  equations  at  each  sample 
point  using  a  Taylor-series  expansion.  This  amounts  to  a  first- 
order  approximation  of  the  required  terms,  with  the  questionable 
assumption  being  that  the  second-  and  higher-order  terms  are  in¬ 
significant.  Additionally,  the  EKF  does  not  accurately  account 
for  the  uncertainty  of  the  underlying  random  variable  in  that  the 
EKF  equations  are  expanded  around  the  a  priori  mean,  with  co- 
variance  expected  to  scale  according  to  the  slope  of  the  function 
at  this  point  only.  The  true  a  posteriori  spread  may  be  signifi¬ 
cantly  different  if  the  function  being  linearized  is  in  fact  quite 
non-linear  in  the  neighborhood  of  the  a  priori  mean.  These  ap¬ 
proximations  may  result  in  large  losses  in  estimation  accuracy 
and  have  been  observed  to  result  in  unstable  filters  [17,18,12]. 


Sigma-point  Kalman  filtering  (SPKF)  is  an  alternate  ap¬ 
proach  to  generalizing  the  Kalman  filter  to  state  estimation  for 
non-linear  systems.  Rather  than  using  Taylor-series  expansions 
to  approximate  the  required  covariance  matrices;  instead,  a  num¬ 
ber  of  function  evaluations  are  performed  whose  results  are  used 
to  compute  an  estimated  covariance  matrix.  This  has  several  ad¬ 
vantages:  (1)  derivatives  do  not  need  to  be  computed  (which 
is  one  of  the  most  error-prone  steps  when  implementing  EKF), 
also  implying  (2)  the  original  functions  do  not  need  to  be  differ¬ 
entiable,  and  (3)  better  covariance  approximations  are  usually 
achieved,  relative  to  EKF,  allowing  for  better  state  estimation, 
(4)  all  with  comparable  computational  complexity  to  EKF. 

SPKF  estimates  the  mean  and  covariance  of  the  output  of  a 
non-linear  function  using  a  small  fixed  number  of  function  eval¬ 
uations.  A  set  of  points  {sigma  points )  is  chosen  to  be  input  to 
the  function  so  that  the  (possibly  weighted)  mean  and  covari¬ 
ance  of  the  points  exactly  matches  the  mean  and  covariance  of 
the  a  priori  random  variable  being  modeled.  These  points  are 
then  passed  through  the  non-linear  function,  resulting  in  a  trans¬ 
formed  set  of  points.  The  a  posteriori  mean  and  covariance  that 
are  sought  are  then  approximated  by  the  mean  and  covariance  of 
these  points.  Note  that  the  sigma  points  comprise  a  fixed  small 
number  of  vectors  that  are  calculated  deterministically — not  like 
the  Monte  Carlo  or  particle  filter  methods. 

Specifically,  if  the  input  random  vector  x  has  dimension  L, 
mean  x,  and  covariance  E*,  then  p  +  1  =  2L  +  1  sigma  points 
are  generated  as  the  set 

X  =  (x,  x  +  X  -  , 

with  columns  of  X  indexed  from  0  to  p ,  and  where  the  matrix 
square  root  R  =  computes  a  result  such  that  E  =  RRT .  Usu¬ 

ally,  the  efficient  Cholesky  decomposition  [19,20]  is  used,  result¬ 
ing  in  lower-triangular  R.  The  reader  can  verify  that  the  weighted 
mean  and  covariance  of  A  equal  the  original  mean  and  covari¬ 
ance  of  random  vector  x  for  a  specific  set  of  {y,  c/C)}  if  we 

define  the  weighted  mean  as  x  =  5Zf=o  the  weighted 

covariance  as  E*  =  J2f=o  ~  —  t)t,  %i  as  the  ith 

column  of  and  both  a and  as  real  scalars  with  the  nec¬ 
essary  (but  not  sufficient)  conditions  that  ^2f=o  =  1  and 

Sf=oa/C)  =  1-  The  various  sigma-point  methods  differ  only 
in  the  choices  taken  for  these  weighting  constants.  Values  for 
the  two  most  common  methods — the  Unscented  Kalman  Filter 
(UKF)  [21-24,18,12]  and  the  Central  Difference  Kalman  Fil¬ 
ter  (CDKF)  [25,26,8] — are  summarized  in  Table  4.  The  UKF 
is  derived  from  the  point  of  view  of  estimating  covariances 
with  data  rather  than  Taylor  series.  The  CDKF  is  derived  quite 
differently — it  uses  Stirling’s  formula  to  approximate  deriva¬ 
tives  rather  than  using  Taylor  series — but  the  final  method  is  es¬ 
sentially  identical.  The  CDKF  has  only  one  “tuning  parameter” 
h,  which  makes  implementation  simpler.  It  also  has  marginally 
higher  theoretic  accuracy  than  UKF  [26],  so  we  focus  on  this 
method  in  the  application  sections  later.  Specifically,  we  used 
h  =  1.2  (which  was  hand-tuned  to  give  somewhat  better  results 
than  the  default  value  of  h  =  >/3),  and  L  =  dim{x}  =  6. 

Before  introducing  the  SPKF  algorithm,  we  reexamine  the 
examples  of  Figs.  2  and  3  using  sigma-point  methods.  In  the 
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Table  4 

Weighting  constants  for  two  sigma-point  methods 


_ y 

UKF  VlTX 
CDKF  h 


a 


(m) 


X 


L+X 

h2-L 


h2 


a 


(m) 


1 


2  (L+X) 
1 


2  h2 


ZTX  +  (1  ~  +  P) 

h2  —  L 
h2 


a 


(c) 


1 


2  (L+X) 
1 


2  h2 


X  =  a2(L  +  k)  —  L  is  a  scaling  parameter,  with  (10-2  <  a  <  1).  Note  that  this 
a  is  different  from  o3m)  and  k  is  either  0  or  3  —  L.  fi  incorporates  prior  infor¬ 
mation.  For  Gaussian  RVs,  ft  =  2.  h  may  take  any  positive  value.  For  Gaussian 
RVs,  h  =  V 3. 


one-dimensional  example  in  Fig.  2,  three  input  sigma  points 
are  needed  and  map  to  the  output  three  sigma  points  shown. 
The  mean  and  variance  of  the  sigma-point  method  is  shown 
as  a  dashed-line  PDF  and  closely  matches  the  true  mean  and 
variance.  For  the  two-dimensional  example  in  Fig.  3,  five  sigma 
points  represent  the  input  random-variable  PDF,  as  shown  in 
frame  (a).  These  five  points  are  transformed  to  the  five  output 
points  in  frame  (b).  We  see  that  the  mean  and  covariance  of  the 
output  sigma  points  (dashed  ellipse)  closely  match  the  true  mean 
and  covariance  and  are  much  better  than  EKF. 

Will  the  sigma-point  method  always  be  so  much  better?  The 
answer  depends  on  the  degree  of  non-linearity  of  the  state  and 
output  equations — the  more  non-linear  the  better  SPKF  should 
be  with  respect  to  EKF.  Given  that  our  battery  model  is  fairly 
linear — the  non-linearities  arise  only  in  the  hysteresis  state  and 
the  OCV  function — we  expect  modest  improvements.  Here,  we 
proceed  to  develop  the  SPKF  algorithm  so  that  we  can  demon¬ 
strate  results  shortly. 

To  use  SPKF  in  an  estimation  problem,  we  first  define  an 
augmented  random  vector  xa  that  combines  the  randomness  of 
the  state,  process  noise,  and  sensor  noise.  This  augmented  vector 
is  used  in  the  estimation  process  as  described  below. 

SPKF  step  1:  State  estimate  time  update.  Each  measure¬ 
ment  interval,  the  state  estimate  time  update  is  computed  by 
first  forming  the  augmented  a  posteriori  state  estimate  vec¬ 
tor  for  the  previous  time  interval:  xk^  =  [(x^)1,  w,  v]T, 

and  the  augmented  a  posteriori  covariance  estimate:  = 

diag(E^_l5  Tiyo ,  T,v).  These  factors  are  used  to  generate  the 
p  ±  1  sigma  points 


X 


a,-\~ 
k- 1 


/v  a,  T 

xk- 1 


+  Y 


a,-\- 
x,k—  1  ’ 


/v  a,  T 

xk- 1  -  Y 


a,  + 
x,k—l 


From  the  augmented  sigma  points,  the  p  +  1  vectors  comprising 
the  state  portion  Xk_x  and  the  p  +  1  vectors  comprising  the 

process-noise  portion  Xkf^  are  extracted.  The  process  equation 

is  evaluated  using  all  pairs  of  Xxk+X  •  and  Xkf^  •  (where  the 
subscript  i  denotes  that  the  ith  vector  is  being  extracted  from  the 
original  set),  yielding  the  a  priori  sigma  points  Xxk~7  for  time 
step  k.  That  is, 

Xi2  =  f(Xi2u,uk-uX?L+u,k-l). 

Finally,  the  a  priori  state  estimate  is  computed  as 


xk  =  E[f(xk-i,  Uk-i,  Wk-i,k  -  l)|Yjt_i] 

p  p 

,k-  1)  = 

1=0  1=0 


E 


E 


2C,+ 

k—  l,z 


Uk- 1 


* k-l,i 


SPKF  step  2:  Error  covariance  time  update.  Using  the  a  pri¬ 
ori  sigma  points  from  step  1 ,  the  a  priori  covariance  estimate  is 
computed  as 


A  Y 

-  X 


k,i 


i) 


T 


SPKF  step  3:  Estimate  system  output  y^.  The  system  output 
is  estimated  by  evaluating  the  model  output  equation  using  the 
sigma  points  describing  the  spread  in  the  state  and  noise  vectors. 
First,  we  compute  the  points  34, i  =  h(Xxkv  ",  Uk,  Xvk'f  t,  k).  The 
output  estimate  is  then 


%  =  E [h{xk,  uk,  vk,  fc)|  Yjt—i] 


p  p 

,  k )  = 

i=0  1=0 

SPKF  step  4:  Estimator  gain  matrix  L^.  To  compute  the  estima¬ 
tor  gain  matrix,  we  must  first  compute  the  required  covariance 
matrices. 


p 

E  «ic)04,i  -  yMj 


1=0 


Then,  we  simply  compute  Lk  =  kT<3 k. 

SPKF  step  6:  Error  covariance  measurement  update.  The 
final  step  is  calculated  directly  from  the  optimal  formulation: 
£+£  =  ^E~k  —  Lk^yfLj.  The  SPKF  solution  is  summarized 
in  Table  5. 


7.  Application  to  battery  management  systems 

7.1.  Cell  and  cell  test  description 

In  order  to  compare  the  various  Kalman  filtering  methods’ 
abilities  to  estimate  SOC  and  SOH,  we  gathered  data  from  a 
prototype  LiPB  cell.  The  cell  comprises  a  LiM^CE  cathode, 
an  artificial  graphite  anode,  is  designed  for  high-power  applica¬ 
tions,  has  a  nominal  capacity  of  7.5  Ah  and  a  nominal  voltage  of 
3.8  V.  SOC  estimation  results  using  EKF  have  been  previously 
reported  for  this  cell,  which  we  call  a  GEN3  cell  [1-6],  and  may 
be  used  as  a  benchmark  for  comparison.  Note  that  the  results  in 
the  companion  paper  [7]  are  for  a  newer  generation  cell,  which 
we  call  a  G4  cell. 

For  the  tests,  we  used  a  Tenney  thermal  chamber  set  at  25  °C 
and  an  Arbin  BT2000  cell  cycler.  Each  channel  of  the  Arbin 
was  capable  of  20  A  current,  and  10  channels  were  connected  in 
parallel  to  achieve  currents  of  up  to  200  A.  The  cycler’s  voltage 
measurement  accuracy  was  ±5  mV  and  its  current  measurement 
accuracy  was  ±200  mA. 

The  cell  test  we  use  here  comprised  a  sequence  of  16  “urban 
dynamometer  driving  schedule”  (UDDS)  cycles,  separated  by 
40  A  discharge  pulses  and  5-min  rests,  and  spread  over  the  90- 
10%  SOC  range.  The  SOC  as  a  function  of  time  is  plotted  in  Fig. 
4(a),  and  rate  as  a  function  of  time  for  one  of  the  UDDS  cycles 
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Table  5 

Summary  of  the  non-linear  sigma-point  Kalman  filter 


Non-linear  state-space  model 

Xk  =  f(xk- 1,  Uk- 1,  Wk-i,k  -  1) 
yk  =  h(xk,  uk,  vk,  k ), 

where  wk  and  vk  are  independent,  Gaussian  noise  processes  of  covariance  matrices  Ew  and  respectively 


Dehnitions:  let 

4  =  WL  4]T’  Xl  =  K**)T.  (^“)T.  (XI)1}1,  p  =  2  X  dime*”) 

Initialization:  for  k  =  0,  set 
vj  =  E[v0] 

£+0  =  E[(*0  -  Sp(*o  -  *o)T] 

Computation:  for  k  =  1,2,...  compute 
State  estimate  time  update 


Error  covariance  time  update 
Output  estimate 

Estimator  gain  matrix 


State  estimate  measurement  update 
Error  covariance  measurement  update 


Vq,+  =  E[xg]  =  [(v+y1,  w,  u]T 
sx!o  =  E[(^o  -  *o+X*o  ~  *o+)T] 

—  o ■> 


'Va'Jr  i  / v<7’+  / pfl,+  i 

Ak-\  —  D4-1*  Xfc-1  +  Yy  Xk-l  Yy^xJ-V 

Xxk-  =f(Xxpil,uk.l,XmkXl,k-l) 

4  =  EL<*T)xL 

Y*  =  EL  aLxL  -  4)(XL  -  v  )T 


yk,,  =  h(xx,uk,xvpu,k) 

-  (mU, 

fk  —  /  j-_n  Otj  34,  i 


%,*  =  EL “!C)(Vu  -  -  »)T 

=  EL^)(-XL  -  w*,-  -  «)r 

=  Vv,/./ Vi 


Xk  =xk  +  Lk(yk  -  yk) 
^t,k  =  ^x,k  -  Lk^y,kLl 


is  plotted  in  Fig.  4(b).  We  see  that  SOC  increases  by  about  5% 
during  each  UDDS  cycle,  but  is  brought  down  about  10%  during 
each  discharge  between  cycles.  The  entire  operating  range  for 
these  cells  (10%  SOC  to  90%  SOC  marked  as  dashed  lines  in 
Fig.  4(a))  is  excited  during  the  cell  test. 

The  data  was  used  to  identify  parameters  of  the  cell  models 
to  be  described  in  the  next  sections.  The  goal  is  to  have  the 
cell  model  output  resemble  the  cell  terminal  voltage  under  load 
as  closely  as  possible,  at  all  times,  when  the  cell  model  input  is 
equal  to  the  cell  current.  Model  fit  was  judged  by  comparing  root- 
mean-squared  (RMS)  estimation  error  (estimation  error  equals 
cell  voltage  minus  model  voltage)  over  the  portions  of  the  cell 
tests  where  SOC  was  between  5%  and  95%.  Model  error  outside 
that  SOC  range  was  not  considered  as  the  HEV  pack  operation 
design  limits  are  10%  SOC  to  90%  SOC. 

7.2.  Enhanced  self  correcting  model  description 

In  order  to  examine  and  compare  performance  of  different 
Kalman  filters,  we  must  first  define  a  discrete-time  state- space 
model  of  the  form  of  (1)  and  (2)  that  applies  to  battery  cells.  Here, 
we  briefly  review  the  “Enhanced  Self-Correcting”  (ESC)  cell 
model  from  references  [3,5].  This  model  includes  effects  due  to 
open-circuit- voltage,  internal  resistance,  voltage  time  constants, 
and  hysteresis. 

State-of-charge  Zk  is  captured  by  one  state  of  the  model.  This 
equation  is 


SOC  trace  versus  time  for  cell  test 


Current  for  one  UDDS  cycle  (zoom) 

80  r1 * i— —  . . . . . . 


—60  * 1 . i  .  .1 . . . . . . . 

205  215  225  235  245 

(b)  Time  (min.) 


Fig.  4.  Plots  showing  SOC  vs.  time  and  rate  vs.  time  for  UDDS  cell  tests.  SOC 
is  shown  in  (a);  rate  for  one  UDDS  cycle  is  shown  in  (b). 
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Empirical  determination  of  OCV  as  a  function  of  SOC 


Empirical  derivative  of  OCV  versus  SOC 


Fig.  5.  Plots  of  (a)  open-circuit- voltage  as  a  function  of  state-of-charge,  and  (b) 
derivative  of  OCV  as  a  function  of  SOC.  In  (b),  raw,  noisy  version  shown  as 
gray,  filtered  derivative  shown  as  black. 


where  rji  is  the  cell  Coulombic  efficiency  at  current  4-1,  A T 
represents  the  inter-sample  period  (in  seconds),  and  C  represents 
the  cell  capacity  (in  As). 

The  time-constants  of  the  cell  voltage  response  are  captured 
by  several  filter  states.  If  we  let  there  be  rif  time  constants,  then 


fk  =  Affk- 1  +  Bfik- 1. 


The  matrix  A  f  e  M"/xw/  may  be  a  diagonal  matrix  with  real¬ 
valued  entries.  If  so,  the  system  is  stable  if  all  entries  have  mag¬ 
nitude  less  than  one.  The  vector  Bf  c  Wxfx  1  may  simply  be  set 
to  rif  “l”s.  The  value  of  rif  and  the  entries  in  the  A  /  matrix  are 
chosen  as  part  of  the  system  identification  procedure  to  best  fit 
the  model  parameters  to  measured  cell  data. 

The  hysteresis  level  is  captured  by  a  single  state 


hk  =  exp 


Qjik-lY  AT 
C 


hk- 1 


+ 


^1  —  exp 


rjiik-\Y  AT 
C 


sgn(/£_i), 


where  y  is  the  hysteresis  rate  constant,  again  found  by  system 
identification. 

The  overall  model  state  is 


*k  =  [ fk  >  hk ,  Zkf- 


The  state  equation  for  the  model  is  formed  by  combining  all  of 
the  individual  equations,  above. 

The  output  equation  that  combines  the  state  values  to  predict 
cell  voltage  is 

yk  =  OCV(zk)  +  Gfk  ~  Rik  +  Mhk, 

where  G  e  Mlx'V  is  a  vector  of  constants  that  blend  the  time- 
constant  states  together  in  the  output,  R  the  cell  resistance  (dif¬ 
ferent  values  may  be  used  for  dis/charge),  and  M  is  the  maximum 
hysteresis  level. 

The  open-circuit- voltage  as  a  function  of  state-of-charge  for 
these  cells  is  plotted  in  Fig.  5(a).  This  is  an  empirical  relationship 
found  by  cell  testing.  For  the  purpose  of  computations  involving 
OCV,  the  final  curve  was  digitized  at  200  points  and  stored  in 
a  table.  Linear  interpolation  is  used  to  look  up  values  in  the 
table. 

The  partial  derivative  of  OCV  with  respect  to  SOC,  required 
by  EKF  but  not  SPKF,  is  plotted  in  Fig.  5(b).  This  relationship 
was  computed  by  first  taking  finite  differences  between  points  in 
the  OCV  plot  in  Fig.  5(a)  and  dividing  by  the  distance  between 
points  (i.e.,  Euler’s  approximation  to  a  derivative).  The  resulting 
data  is  too  noisy  to  be  of  practical  use,  as  shown  in  the  gray 
line  of  Fig.  5(b).  It  was  filtered  using  a  zero-phase  low-pass 
filter,  resulting  in  the  black  line  of  Fig.  5(b),  which  is  used  in 
EKF.  This  relationship  is  also  digitized  at  200  points,  and  linear 

Estimated  and  actual  cell  terminal  voltage 
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Fig.  6.  Voltage  prediction  using  cell  model.  Gray  is  true  voltage,  black  is  esti¬ 
mated  voltage.  In  (b)  a  zoom  of  voltage  prediction  for  one  UDDS  cycle  at  around 
50%  SOC  is  shown. 
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interpolation  into  the  table  of  values  is  used  when  computations 
requiring  this  function  are  performed. 

Other  parameters  are  fit  to  the  cell  model  using  a  method 
previously  described  ([5],  Section  4).  In  particular,  the  model 
employs  four  low-pass  filter  states  (rif  =  4),  a  nominal  ca¬ 
pacity  of  7.5  Ah,  and  an  inter-sample  interval  of  AT  =  1  s. 
Further,  Af  =  diag{0.9624,  0.8509,  0.9981,  0.99999},  Bf  = 
1CT4  x  [1,  1,  1,  1]T  (chosen  to  scale  all  states — including  the 
hysteresis  and  SOC  state — to  roughly  the  same  dynamic  range), 
G  =  [-0.5256,  -1.3258,  -0.1855,  0.0012],  y  =  2.2523,  rj  = 
1,  M  =  74.7  mV,  the  charging  resistance  was  2.6  m^,  and  the 
discharging  resistance  was  2.7  mQ.  These  are  the  same  values 
used  to  generate  the  results  in  [5].  There  is  very  close  agree¬ 
ment  between  the  cell  model  voltage  prediction  and  the  cell  true 
voltage.  This  is  illustrated  in  Fig.  6(a).  To  better  illustrate  the 
model’s  fidelity,  refer  to  the  zoom  on  one  UDDS  cycle  in  the 
50%  SOC  region,  shown  in  Fig.  6(b). 

7.3.  Examples  of  EKE  versus  SPKF 

Using  the  data  fit  to  the  ESC  model,  we  ran  an  EKF  and  an 
SPKF  to  predict  SOC.  Note  that  the  EKF  results  in  this  section 
are  identical  to  those  in  [6]  and  are  replicated  here  for  clarity  of 
comparison  with  the  SPKF  results.  All  “tuning”  parameters  were 
the  same  between  EKF  and  SPKF;  the  model  and  all  cell-test 
data  files  used  were  identical. 


Estimation  error:  ESC,  nf=  4 


Estimation  error:  SPKF  ESC,  nf- 4 


Fig.  7.  SOC  estimation  error  for:  (a)  EKF  vs.  (b)  SPKF,  when  correctly  initial¬ 
ized. 


Estimation  error:  ESC,  nf= 4 


Estimation  error:  SPKF  ESC,  rif=4 


(b)  Time  (min.) 


Fig.  8.  SOC  estimation  error  for:  (a)  EKF  vs.  (b)  SPKF,  when  not  correctly 
initialized. 

SOC  estimation  error  is  plotted  in  Fig.  7  for  both  technologies 
when  the  filters  are  correctly  initialized  (to  initial  SOC  =100%). 
The  EKF  RMS  SOC  estimation  error  is  0.64%  while  the  SPKF 
RMS  SOC  estimation  error  is  0.49% — an  improvement  of  over 
23%.  The  maximum  absolute  SOC  estimation  error  was  1.10% 
for  EKF  and  0.90%  for  SPKF — an  improvement  of  18%.  The 
bounds  correctly  encompassed  the  zero  point  95. 1 1  %  of  the  time 
for  SPKF,  but  only  74.73%  of  the  time  for  EKF. 

In  Fig.  8,  we  compare  the  results  for  EKF  versus  SPKF  when 
the  filters  were  incorrectly  initialized  to  SOC=80%  rather  than 
the  correct  100%.  Both  filters  had  equivalent  uncertainty  matri¬ 
ces  on  the  states,  so  recovery  from  the  initial  error  took  a  similar 
amount  of  time.  Error  bounds  for  the  SPKF  method  were  slightly 
tighter,  however. 

The  EKF  had  an  RMS  SOC  estimation  error  of  0.75%,  while 
the  SPKF  had  an  RMS  SOC  estimation  error  of  0.69%.  This  is  an 
8%  improvement,  which  at  first  appears  small  compared  to  the 
first  trial  case.  Notice  that  the  RMS  error  is  here  dominated  by  the 
initial  convergence  of  the  estimators,  which  is  nearly  equivalent 
for  both  technologies.  The  lower  error  of  SPKF  may  be  attributed 
mostly  to  its  convergence  closer  to  the  true  SOC  after  the  initial 
transient.  Also,  the  SPKF  bounds  were  correct  97.86%  of  the 
time,  while  the  EKF  bounds  were  correct  only  93.86%  of  the 
time.  These  numeric  results  are  replicated  in  Table  6  for  easy 
comparison.  In  steady  state  we  observe  SPKF  to  improve  RMS 
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Table  6 

Comparison  of  EKF  vs.  SPKF  in  UDDS  test  results  predicting  SOC 
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Correctly  initialized 

Incorrectly  initialized 

RMS  error  (%) 

Maximum  error  (%) 

Bounds  error  (%) 

RMS  error  (%) 

Bounds  error  (%) 

EKF 

0.64 

1.10 

25.27 

0.75 

6.14 

SPKF 

0.49 

0.90 

4.89 

0.69 

2.14 

Improvement 

23 

18 

81 

8 

65 

error  by  about  20%  and  to  greatly  improve  prediction  of  error 
bounds  because  of  a  better  state  error  covariance  estimate. 


We  first  derive  a  fundamental  equation  governing  Gaussian 
estimation,  and  then  apply  it  to  the  problem  at  hand. 


8.  Conclusions 

This  paper  has  considered  several  methods  for  state  estima¬ 
tion  of  a  battery  cell  with  application  to  battery-management  sys¬ 
tems  of  hybrid- electric  vehicles.  These  algorithms  are  based  on 
optimal  estimation  theory  (also  known  as  sequential  probabilis¬ 
tic  inference)  and  encompass  several  members  of  the  Kalman 
filter  family.  In  particular,  we  have  shown  how  the  Kalman  filter, 
the  extended  Kalman  filter  and  the  sigma-point  Kalman  filters 
can  be  derived  and  applied  to  state  estimation. 

In  prior  work  we  have  shown  how  the  extended  Kalman  fil¬ 
ter  can  be  used  in  battery-management  systems.  Here,  we  have 
shown  that  the  sigma-point  Kalman  filter  may  also  be  used  for 
LiPB-cell  state  estimation,  at  an  equivalent  computational  com¬ 
plexity.  Data  from  testing  using  a  third-generation  prototype 
LiPB  cell  shows  that  SPKF  gives  superior  results.  The  RMS 
SOC  estimation  error  is  lower,  the  maximum  SOC  estimation 
error  is  lower,  and  the  error  bounds  produced  on  the  estimate  are 
more  accurate. 

The  state  estimate  produced  by  SPKF  can  also  be  used  to 
predict  available  power  [27],  or  to  effect  equalization  via  SOC 
[6] .  Furthermore,  the  estimate  can  be  made  more  accurate  if  the 
cell  parameters  are  also  estimated  in  real  time  to  account  for  any 
manufacturing  differences  between  cells,  and  to  track  the  effects 
of  aging.  Details  on  how  to  do  this  using  SPKF  are  discussed  in 
the  companion  to  this  paper  [7] . 


Theorem.  If  x  and  y  are  jointly  Gaussian  vectors  with  means 
x  and  y  and  joint  covariance  E,  then  p(x\y)  is  Gaussian  with 
mean  x  +  E^-yE”1  (y  —  y)  and  covariance  Exx  —  E^E^1  Eyx, 

and  thus  E[v|y]  =  x  +  Hxy  E”,1  (y  —  y). 


Proof.  First,  we  write, 


/  i  x  p(x>y) 
p(x  |  y)  =  — — —  a  exp 

piy) 
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where  E  = 


^xx  E 
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is  the  joint  covariance  matrix.  Then, 


substitute  the  transformation 
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maintaining  the  inner  matrix  without  expansion.  The  terms  in 
the  exponent  of  (12)  become 
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Appendix  A.  Proof  of  Gaussian  recursion 

In  this  appendix,  we  prove  the  claims  of  Eqs.  (3)  and  (4) 
using  the  definitions  of  Eqs.  (5)— (1 1).  That  is,  we  find  E[jCfc|Yfc] 
under  the  basic  assumption  that  all  densities  remain  Gaussian. 


verifying  the  claim. 

We  can  now  proceed  to  verify  Eqs.  (3)  and  (4)  by  find¬ 
ing  E[jc*|Y*].  First,  we  define  xf  =  Xk  —  xf  where  xf  = 
E[*jt|Y*_i]  and  %  =  yk  -  yk  where  yk  =  E[^|Y^_!].  By  as- 
sumption,  xf  and  y^  are  jointly  Gaussian,  and  by  construction 
have  zero  mean.  Then,  by  Theorem  1  we  have  that  E[xf  |y&]  = 

^xy,k^fk%,  which  we  define  to  be  Lk%. 

Secondly,  we  must  show  that  E[x^|Y^_i]  =  E[ji^~]  and 
therefore  xjj  is  uncorrelated  with  Y^_i  (independent  because 
Gaussian).  This  is  straightforward  by  substitution: 

E[v^|Y^_i]  =  E[xk  -  E{v^|Y^_i}|Y^_i]  =  0  =  E[xf]. 
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Therefore,  we  can  write 

E[5^|Y*]  =  E[xk  |  Y*_i,  yk] 

=  E[xk \yk]  =  E[xk\yk]  -  E[x^\ %] . 

Lk(yk-yk)  x+ 

From  this  last  line,  we  solve  for  the  a  posteriori  state  estimate 
V  =  V  +  Lk(yk  -  %),  which  verifies  Eq.  (3)  □. 

The  covariance  of  xk  may  be  computed  using  Eq.  (3). 

£ f  =  E[{(x&  -  xk)  -  Lkyk}{(xk  -  xk)  -  Lkyk}T] 

=  ~  Lk^[yk(xk)  ]Lk  ~  %  ]  Lk  +  LkYuLk 

~  ^x,k  ~  Lk'EytkLj, 

which  verifies  Eq.  (4). 
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